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Abstract — Block tridiagonal systems play a key role in all 
Kalman smoothing applications, including the classic Rauch- 
Tung-Striebel smoother, as well as more modern variants that 
incorporate nonlinear models, inequality constraints on the state, 
and robust penalties on state and measurement components. 

The paper begins with a fundamental result that connects 
the condition number of a common block tridiagonal system 
to properties of the individual blocks. As a consequence, we 
obtain sufficient conditions for the stability of Kalman smoothing 
formulations. The result also illustrates how unstable Kalman 
smoothing problems can arise, and how to easily stabilize them. 

Then, turning our attention to algorithms, we show that the 
classic Rauch-Tung-Striebel smoother is an implementation of 
the forward Thomas algorithm for block tridiagonal systems. 
We reveal a flaw in the existing theory for characterizing the 
numerical stability of the algorithm, and provide a new numerical 
stability result that shows the condition number of every recursively 
modified block is bounded by the condition number of the whole 
system as the algorithm proceeds. 

Finally, we present a new backward block tridiagonal Thomas 
algorithm, and show that the condition numbers of the recursively 
modified blocks in this algorithm are independent of the condition 
number of the full block tridiagonal system. 



I. Introduction 

Kalman filtering and smoothing methods form a broad 
category of computational algorithms used for inference on 
noisy dynamical systems. Since their invention 031 and early 
development fT6lL these algorithms have become a gold stan- 
dard in a range of applications, including space exploration, 
missile guidance systems, general tracking and navigation, 
and weather prediction. In 2009, Rudolf Kalman received the 
National Medal of Science from President Obama for the 
invention of the Kalman filter. Numerous books and papers 
have been written on these methods and their extensions, 
addressing modifications for use in nonlinear systems lfT4l . 
robustification against unknown models l22i smoothing data 
over time intervals Q, ((J, Kalman smoothing with unknown 
parameters |8], constrained Kalman smoothing (7), and im- 
proving robustness to bad measurements fTTlL EH, fTOlL l24ll 
|20L flgl. flTl. 151 and sudden changes in state l3l HI. fl3l 
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Therefore, we use the term 'Kalman smoothing' to refer to 
any system that can be represented graphically by Figure [T] 
In order to motivate the analysis in this paper, we begin by 
reviewing the classic Kalman smoothing problem, and show 
that it reduces to solving a linear system with block tridiagonal 
structure. 

The model corresponding to Figure [T] is specified as follows: 

X\ = gi(>o)+wi, 

Xk = gk(xk-i)+Wk k = 2,...,N, (1.1) 
Zk = h k (x k )+v k k=l,...,N, 

where Wk, v k are mutually independent random variables 
with known positive definite co variance matrices Q k and R k , 
respectively. Extensions have been proposed that work with 
singular Q k and R k (see e.g. fT9lL l9l), but in this paper we 
confine our attention to the classic nonsingular case. 

We have x k ,w k G W 1 , and z k ,v k £ W 1 ^ , so measurement 
dimensions can vary between time points. The classic case is 
obtained by making the following assumptions: 

1) xq is known, and g k , h k are known linear functions, 
which we denote by 

8k (xk- 1 ) = GkXk- i h (x k ) = H k x k (1.2) 

where G k G R nxn and H k G R m ( k ) xn , 

2) w k , v k are mutually independent Gaussian random vari- 
ables. 

If the models are linear and the noises are Gaussian, using 
B ayes' theorem, we have 

P {x k \zk) « P (zk\xk)p(x k ) = p(v k )p(w k ) , (1.3) 

and therefore the likelihood of the entire state sequence {x k } 
given the entire measurement sequence {z k } is proportional to 

N N 1 

n p(^)p( w ^) n ex p ( ~ o (** ~ H k*k) T R k 1 (z k - H k x k ) 

k=l k=l Z 

- ~(x k - G k x k _i) T Q~ l (x k - G k x k _i)^ . 

(1.4) 

A better (equivalent) formulation to ( |I.4| ) is minimizing its 
negative log posterior: 

N 1 

min/({**}) := £ - (zk - H k x k ) T R^ 1 (z k -H k x k ) 

w k=\ z (1.5) 

+ - (** - G k Xk-\ ) T Q k l (xk ~ G k Xk-i) • 
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Fig. 1. 



Dynamic systems amenable to Kalman smoothing methods. 



To simplify the problem, we now introduce data structures 
that capture the entire state sequence, measurement sequence, 
covariance matrices, and initial conditions. 

Given a sequence of column vectors {u k } and matrices {7^} 
we use the notation 

~7i ••• 



vec({v^}) 



, diag({7i}) 













T 2 





T N 



We make the following definitions: 

R = diag({/?jk}) x = vec ({**}) 
Q = dmg({Q k }) C = vec({x ,0,. 
H = dmg({H k }) z = vec({zi,Z2, 
I 



,0}) 



-G 2 I 



-Gn 



(1.6) 



(1.7) 



With definitions in 
written 



L6| ) and ( |I77| ), problem ( [13] ) can be 



mmf(x) 



\Hx 



-z\\U 



\q-i > 



(1.8) 



where \\a\\]^ = a 1 Ma. It is well known that the MAP was a 
least-squares problem, but this derivation makes the structure 
fully transparent. We now write down the linear system that 
needs to be solved in order to find the solution to this least- 
squares problem: 



(H T R- l H + G T Q- l G)x 



H T R~ l z 



G T Q-^. (1.9) 



The linear system in ( |L9| ) has a very special structure: it 
is symmetric positive definite block tridiagonal, since G is 
nonsingular and Q is positive definite. Direct computation 
shows that 

|Z>i 



C 



(H T R~ 



G T Q~ 



A 2 




D 2 A 







A N D A 



(1.10) 



with A k G . 
A k ~- 



and D k G \ 



defined as follows: 



G M Q M G k +i 



H k R l 



G k 



(I.H) 



This block tridiagonal structure was noted early on in l25ll 
lfT2l . l26l . These systems also arise in many recent extended 
formulations, see e.g. (15)], Q (12)], E (8.13)]. 

The paper proceeds as follows. Section [TT] focuses on the 
subclass of block tridiagonal systems that capture all of the 
systems of interest to various Kalman smoothing applications. 
Bounds on the eigenvalues of these systems are obtained in 
terms of the behavior of the individual components. It is shown 
how these bounds are related to the numerical stability of 



Kalman smoothing formulations. In Section III it is shown 
that the classic RTS smoother implements the forward block 
tridiagonal Thomas algorithm for the system ( |L9| ). We demon- 
strate a flaw in the stability analysis for this algorithm given 
in 0, and prove a new more powerful stability result that 
shows any stable system from Section [TT] can be solved using 
the forward block Thomas algorithm. Finally, in Section [V| we 
introduce a new backward Thomas algorithm, and show that it 
behaves stably independent of the condition number of the full 
system. We conclude with a discussion of the consequences 
of the new results. 

II. Characterizing block tridiagonal systems 
Consider systems of form 

B = a T qa (II. 1) 

where 

^10 



^ = diag{^i,...^}, 



a = 



a>2 



I 



Let Amin, Amax, and <7 m in, <7 max denote the minimum and max- 
imum eigenvalues and singular values, respectively. Simple 
upper bounds on the lower and upper eigenvalues of B are 
easily obtained. 



Theorem 2.1: 

^min(^)^min(^) < ^min(^) < Ksaiifi) < A max (^)(T^ ax (a) . 

(H.2) 

Proof: For the upper bound, note that for any vector v, 

v T a T qav < A max (^) \\av\\ 2 < ^max(^)^ ax (^) ||v|| 2 . 

Applying this inequality to a unit eigenvector for the maximum 
eigenvalue of c gives the result. The lower bound is obtained 
analogously: 

v T a T qav > K^{q) \\av\\ 2 > (a) ||v|| 2 . 

Applying this inequality to a unit eigenvector for the minimum 
eigenvalue of c completes the proof. ■ 
From this theorem, we get a simple bound on the condition 
number of k(B): 

k ^ = ^max(^) < j^ax(g)g^ax(g) ^ 

Since we typically have bounds on the eigenvalues of q, all 
that remains is to characterize the singular values of a in terms 
of the individual a k . This is done in the next result which uses 
the relation 



(i- 



T 

a a - 



-a\a 2 



CL2 



-ajai 







(H.4) 



where we define a^+i := 0, so that the bottom right entry is 
the identity matrix. 

Theorem 2.2: The following bounds hold for the singular 
values of a: 

^2 



min { 1 

k 1 



< 



ol^{a T a 



max { 1 

k 1 



(%+i)- 
< 

(0*+iH 



-'max 



-'max 



0^+1 )} 



Smax(^) < 
Omax(^) + Omax(0*+l)} 



(115) 



Proof: Let v = vec({vi, . . . , v^}) be any eigenvector of 
a T a, so that 

Av . (II.6) 



T 

a av - 



Without loss of generality, suppose that the v k component 
has largest norm out of [1,...,A/]. Then from the kth block 
of jll.6| ), we get 



ajcVk-i + ( / + ^+i^m) v £ + ^+i v £+i = Xv k . (117) 
Let u k = Multiplying \IIJ\ on the left by vj, dividing by 
||v^|| 2 , and rearranging terms, we get 



1 



- Ukd k+x Clk+\Uk - 



-u k a k - 



T v k+l 

u k a k ^ - 



Wl " K+1 ||v.|| 

S ^max 

This relationships in (|II.8|) yield the upper bound 



A < 1 + CT max (^ + i) + Omax(^) + ^max(^+l) 



and the lower bound 

A > 1 + u k aj 

> 1 + o£in(a*+l) " 



U k a k +ia k +\U k — c max 
^2 



- cr max (^ + i) 
(fit) - a ma x(a<; + i) . 



(H.8) 



(H.9) 



(11.10) 



Taking the minimum over k in the lower bound and maximum 
over k for the upper bound completes the proof. ■ 
Corollary 2.1: Let v mm be the eigenvector corresponding to 
A rmn (^ T a), and suppose that ||v^ m || is the component with the 
largest norm. Then we have the lower bound 

Amin^V) > l + a min (^ + i)-CJ max (^)-0- max (^ + i) . (11.11) 

In particular, since a^^i = 0, 

A min (a T a) > 1 - o max (a N ) . 



(11.12) 



The bound 



11.12 



reveals the vulnerability of the system a T a to 
the behavior of the last component. 

For Kalman smoothing problems, the matrix a T qa corre- 
sponds to G T Q~ l G. The components a k correspond to the 
process models G k . These components are often identical for 
all k, or they are all constructed from ODE discretizations, 
so that their singular values are similarly behaved across k. 
By ( |II.11| ), we see that the last component emerges as the 
weakest link, since regardless of how well the G k are behaved 
for k = 2, . . . , N — 1, the condition number can go to infinity 
if any singular value for Gn is larger than 1. The last block 
can thus have a very powerful effect on the entire system. 

III. Forward Block Thomas Algorithm 

We now present the forward block Thomas algorithm. 
Forward here refers to the fact that the algorithm starts and 
ends in the top left corner, in other words, it goes forward 
and then backward. In contrast, the backward block Thomas 
algorithm we present in Section |V| will start and end in the 
bottom right corner, so it goes backward, then forward. 

Suppose for k = 1, . . . ,N, b k £ R nxn , e k £ R nxi , r k £ B nxi , 
and for k = 2, . . . ,N, c k £ R nxn . We define the corresponding 
block tridiagonal system of equations 



(b x c\ 
c 2 b 2 





Vo 







cn-i 




o\ / e x \ ( n \ 



b N _i 
c N 



b N J 



eN-l 
e N 



rN-l 
\ rN ) 



(IIL1) 



For positive definite systems, the forward block Thomas 
algorithm is defined as follows |8, algorithm 4]: 
Algorithm 3.1: The inputs to this algorithm are {c k } k=2 , 

{M*=i' and Mf=i where each c k € RnXn > b k ^ R nxn , and 
r k £R nxi . The outputs are a sequence {e k }^ =l and 7] where 
each e k £ n nxi and 7] £ R. The sequence solves equation ( III. 1 



and the A that is equal to the log of the determinant of the 



block tridiagonal matrix in equation ( III. 1 

1) Set d\ = b\. Set si = r\. 

2) For k = 2 To N : 

• Set d k = b k 

• Set s k = r k 

3) Set e^ = d^ l SN. 

4) For k = N- \ To 1 

• Set e k = d^ l (s k 

5) Set 7] 



-c k d k \c\. 
-c k d k \s k -i. 



:Lt!logdet(4). 



Before we discuss stability results for this algorithm (see 



theorem 3.2), we prove that the RTS smoother is an imple- 
mentation of this algorithm for matrix C in ( |I.10| ). 

Theorem 3.1: When applied to C in ( |I10| ), Algorithm 3.1 
is equivalent to the Rauch-Tung-Striebel (23ll smoother. 

Proof: Looking at the very first block, we substitute in 
the Kalman data structures into step 2 of Algorithm [37T 



Writing this step requires introducing some structures which 
may be familiar to the reader from Kalman filtering literature. 

p^-^q^+hJr^ 1 ^ 



1 21 1 



PT 1 - 

r 2 2 ' 



(GiP 1 . 1 Gf + G2)" 



P^+GjQ~ l G 2 y (q~ 1 G 2 ) (IIL2) 



r 2|l 



d 2 = c 2 — a^d l l a 2 = P r 1 



2 |2 TU 3 ^3 G 3 



These relationships can be seen quickly from (21 Theorem 
2.2.7]. The matrices P k \ k , Pk\k-\ often appear in Kalman 
literature: they represent covariances of the state at time k 
given the the measurements and the covariance 

of the a priori state estimate at time k given measurements 
{zi , • • • , Zk-i } , respec tively . 
The key fact from |III.2| is that 



d 2 



: P 2\2 



gJq^g 3 . 



Using the same computation for the generic tuple 
rather than (1,2) establishes 



d k = P l 



k\k 



(III.3) 



We now apply this technique to the right hand side r 
H T R- l z + G T Q- l w. We have 

:= (Q 2 l G 2 ) T (p^+GlQ- l G 2 y l {hJr-'z^gJp-^ 



yi\\ 
yi\i 



n 2 R 2 l z 2 +y 2 \i 



s 2 = r 2 — ajd l 1 r\ 



'-yi\i 



(III.4) 

These relationships also follow from O Theorem 2.2.7]. The 
quantities j 2 |l an d J2|2 ma Y be familiar to the reader from 
the information filtering literature: they are preconditioned 
estimates 

yk\k Pjc\k x k\k > yk\k-i P k \k-\ x k\k-i > (iii.5) 

where x k \ k is the estimate of %k given {zi, . . . and 

x k\k-\ — GkXk-\\k-\ 

is the best prediction of the state Xk given {z\ : . . . 

Again, by applying the computation to we have precisely 
that S£ = y k \ k . From these results, it immediately follows that 
eN computed in step 3 of Algorithm 3J_ is the Kalman filter 
estimate (and the RTS smoother estimate) for time point N 
(see gjO) ): 

-l 



'N\N X N\N — X N\N ■ 



(III.6) 



We turn now to Step 4 of Algorithm 3.1 First, follow- 
ing (23] (3.29)], we define 



To save space, we also use shorthand 



P k :=Pi 



k\k, Xk'-—*k\k- 



(HI.7) 



(III.8) 



At the first step, we obtain 



eN-i 



(Pn-i 



(s N -\ -c T N e N ) 



-GjjQ N l G N ) l (P N !_ l x N -i-GjjQ N l x N ) 



G^Q n 1 Gn) 



1 P n _iXn-i 



= X N -\ —CN-i(G n x N -i —x N ) 

= x iV- 1 L/V- 1 + C/V- 1 ( X N\N ~ GnX n _ , 

(HI.9) 

where the Woodbury inversion formula was used to get from 
line 3 to line 4. Comparing this to (23l (3.28)], we find 
that eN-i = x N _i\ N , i.e. the Rauch-Tung-Striebel smoothed 
estimate. The computations above, when applied to the general 
tuple (&,/:+ 1) instead of (N—1,N), show that every is 
equivalent to x k \ N , which completes the proof. ■ 

In 1965, Rauch, Tung and Striebel showed that their 
smoother solves the maximum likelihood problem for 
P({*fc}|fe}) E3' which is equivalent to ( |T3] ). Theorem 3.1 
adds to this understanding, showing that it is precisely the for- 
ward block Thomas algorithm. Moreover, the proof explicitly 
shows how the Kalman filter estimates x k \ k can be obtained as 
the forward block Thomas proceeds to solve system ( |L9| ). For 
efficiency, we did not include any expressions related to the 
Kalman gain; the interested reader can find these relationships 
in Chapter 2]. They emerge through application of the 
Woodbury matrix identity to the equations we presented. 

We now turn out attention to the stability of the forward 
Thomas algorithm for block tridiagonal systems. One such 
analysis appears in [ 8 , Lemma 6], and has been used frequently 
to justify the forward Thomas as the method of choice in 
many Kalman smoothing applications. Below, we review this 
result, and show that it has a critical flaw precisely for Kalman 
smoothing systems. 

Lemma 3.1: Suppose we are given sequences {qk}k=o 9 
{ak}k=2> anc * i u k}k=v wn ere each q k G R nxn is symmetric 
postive definite, each Uk G R nxn is symmetric semi-postive 
definite, each a k G R nxn . Define b k G R nxn by 



u k + q k 



» wnere k 



l,...,iV 



Define c k G R nxn by 

Ck = q k l dk ^ where k = 2,...,N 

Suppose there is an a > such that all the eigenvalues of 
ajq^ l ak are greater than or equal a. Suppose there is a j3 > 
such that all the eigenvalues of bk are less than or equal /3. 



Further suppose we execute Algorithm 3.1 with corresponding 



input sequences {bk} k=l and {ck} k=2 - It follows that each dk 
generated by the algorithm is symmetric positive definite and 
has condition number less than or equal f3/a. 



To understand what can go wrong with this algorithm, 
consider in the Kalman smoothing context where the 
corresponding matrix entries are given by | 8, equations (12) 
and (13)]. The matrix in Lemma [TT] is given by 



and 



T -1 

-a N j rl q NJrl a^j r \ 



where the correspondence to (8j (13)] is given by b^ corre- 
sponds to Bn, corresponds to Qn, corresponds to Gn, 
and un corresponds to 

In the context of |8], G k +\ is the model for the next state vector 
and at k = N there is no next state vector. Hence, G^+i = 0. 
Thus, in the context of Lemma |3J"] = and hence a = 
which contradicts the Lemma assumptions. 

Does this mean that the forward block Thomas algorithm 
(and hence the Kalman filter and RTS smoother) are inherently 
numerically unstable, unless they have measurements to stabi- 
lize them? It turns out that this is not the case; in fact we can 
prove a powerful theorem that relates numerical stability of 
the forward block Thomas algorithm to the condition number 
of the system ( |II.1| ), already characterized in Theorem 2.1 



Theorem 3.2: Consider any block tridiagonal system A G 
R N of form ( |III.1| ). and suppose we are given a lower bound 
(Xl and an upper bound ajj on the eigenvalues of this system: 



If we apply the block Thomas iteration 



(III. 10) 



then 



d k = b k -c k d k \cj, 



0<a L < KtmiAk) < ^(4) < V£ . 



(III.ll) 



In other words, the block Thomas iteration preserves eigen- 
value bounds (and hence the condition number) for each block, 
and hence will be numerically stable when the full system is 
well conditioned. 

Proof: For simplicity, we will focus only on the lower 
bound, since the same arguments apply for the upper bound. 
Note that b\ =d\, and the eigenvalues of d\ must satisfy 

OCl < ^min(^l) 

since otherwise we can produce an unit-norm eigenvector vi G 
W 1 of d\ with v\d\v\ < CCl, and then form the augmented unit 
vector v\ G M, N with vi in the first block, and every other entry 
0. Then we have 

vjAvi < aL , 



which violates ( |V.2| ). Next, note that 
fbi 

d 2 



S ] AS 





Vo 



CN-l 





where 



b N _i 
c N 

1 T 



(III. 12) 



-c 2 d: 



Si = 





/ 



V 








Suppose now that d 2 has an eigenvalue that is less than a^. 
Then we can produce a unit eigenvector v 2 of d 2 with vjd 2 v 2 < 
(Xl, and create an augmented unit vector 

V2 = [0\ xn v\ lxn ( N _ 2 )] T 

which satisfies 

vlSiAS[v 2 < a L . (III. 13) 

Next, note that 

v| := vJSi = [-vjc 2 d~ l vj lxn{N _ 2) ] T , 
so in particular ||v2|| > 1. From ( |III.13| ), we now have 

v 2 Av 2 <a L < a L \\v 2 \\ , 

which violates ( |V.2| ). To complete the proof, note that the 
lower n(N— l)xn(N— 1) block of S\ASj has the same form as 



A, with ( |V.2| ) holding for this modified system. The reduction 
technique can now be repeatedly applied. ■ 
Note that Theorem |3.2| applies to any block tridiagonal 
system. When applied to the Kalman smoothing setting, if 
the system G T Q~ l G is well-conditioned, we know that the 
forward block Thomas (and hence the Kalman filter and 
RTS smoother) will behave well for any measurement model. 
Moreover, even if G T Q~ l G has a bad condition number, it is 
possible that the measurement term H T R~ l H (see ( |I.9| Q will 
improve the condition number. More general Kalman smooth- 
ing applications may not have this advantage. For example, the 
initialization procedure in [ 5] requires the inversion of systems 
analogous to G T Q~ l G. 

IV. INVERTIBLE MEASUREMENT COMPONENT 



We now return to the system ( |I.9| ), and consider the case 
where H T R~ l H is an invertible matrix. Note that this is not 
true in general, and in fact our analysis, as applied to the 
Kalman smoothing problem, did not use any assumptions on 
this term. However, if we know that 



F \=H T R~ l H 



(IV.l) 



is an invertible matrix, then we can consider an alternative 
approach to solving ( |L9| ). Applying the Woodbury inversion 
formula, we obtain 

(G T Q~ l G + F)~ l =F~ l -F- l G T (Q + GF- l G T )- l GF~ l 

(IV.2) 

Now, the solution to ( |L9| ) can be found by applying this 
explicit inverse to the right hand side 

rhs =H T R- l z + G T Qr l t; 



d 2 = b 2 - c 2 d x c 2 



and the key computation becomes 

(2 + GF~ l G T )x = GF~ l rhs . 



(IV.3) 



Note that the matrix Q + GF~ l G T is block tridiagonal, since 
Q,F~ l are block diagonal, and G is lower block bidiagonal. 
Therefore, we have reduced the problem to a system of the 
form Moreover, at a glance we can see the lower eigen- 

values of this system are bounded below by the eigenvalues of 
<2, while upper bounds can be constructed from eigenvalues of 
2, Gk and F~ l . Under very mild conditions, this system can 
be solved in a stable way by the forward tridiagonal algorithm, 
which would also give a modified filter and smoother. This is 
not surprising, since we have assumed the extra hypothesis 
that F is invertible. 

V. Backward Block Tridiagonal Algorithm 

Now that we have abstracted Kalman smoothing problems 
and algorithms to solutions of block tridiagonal systems, we 
can consider novel techniques for their solution. In this section, 
we propose a backward block Thomas algorithm, and then 
show it has a very interesting property when applied to the 
Kalman smoothing setting. Recall that backward means the 
algorithm starts and ends in the lower right hand corner, so it 
first goes backward, and then forward. 

To derive the backward algorithm, we begin with system 
( |III.1| ). Subtracting c N b^ 1 times row N from row N—l, we 
obtain the following equivalent equation: 



(bx 


T 











b 2 
























bN-2 


c N-l 










cn-i 


bn-l -cJfb^ l c N 





Vo 







C N 


bj 



\ e N / 



rN-l 



n 



-cJjb N l r N 



J 



\ r N 

We iterate this procedure until we reach the first row of the 
matrix, using d k to denote the resulting diagonal blocks, and 
s k the corresponding right hand side of the equations; i.e., 

d N = b N , d k = b k - cJ+id^Ck+i (k = N- 1, • • • , 1) 

s N = e N , s k = r k - (* = N ~ 1 j ' ' ' j !) • 

We obtain the following lower triangluar system: 

0\ 



(di 

C2 





d 2 







cn-i 







d N -\ 
\ c N b N J 

Now we can solve for the first block vector and then proceed 
back down, doing back substitution. We thus obtain the 
following algorithm: 

Algorithm 5.1: A block tridiagonal extension of the back- 
ward Thomas algorithm: The inputs to this algorithm are 
{c k }, {b k }, and {r k }. The output is a sequence {e k } that 
solves equation the 7] which is equal to the log of the 

determinant of the block tridiagonal matrix in equation 



( 61 \ 



\ e N J 



( 51 \ 

S2 



S N -l 



(V.l) 



1) Set dN = btf and sn = r#. 

2) For k = N — 1 , . . . , 1, set a\ = bk — cj^ d~ l 



r k- c k+\ d k+l S k+l' 

3) Set e x =d~ l s x . 



4) For k = 2,...,N, set e k = d k l (s k -c k e k -i). 

5) Set T]= E^^og det(4). 

First, we show that the backward block Thomas algorithm 
has the same numerical stability result as the forward block 
Thomas algorithm. 

Theorem 5.1: Consider any block tridiagonal system A G 
R N of form and suppose we are given the bounds 

and Oiu for the lower and upper bounds of the eigenvalues of 
this system 

< a L < A m i n (A) < Amax(A) < ecu . (V.2) 
If we apply the backward block Thomas iteration 

dk = h 
then we have 

0<a L <A min (4)<A max (4)<a t / Vfc. (V.3) 
Proof: Note first that d^ = b^, and satisfies (|V.3|) by the 



" c J+l^+l c A:+l 



same argument as in the proof of Theorem |5.1| Define 




and note that 
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Now an analogous proof to that of Theorem [571 
applied to show the upper n(N—l)xn(N—l 
satisfies ( |V.3| ). Applying this reduction iteratively completes 
the proof. ■ 

Theorems 13.21 and I57T1 show that both forward and backward 
Thomas algorithms are numerically stable when the block 
tridiagonal systems they are applied to are well conditioned. 
However, a different analysis can also be done for a particular 
class of block tridiagonal systems. This result, which applies 
to Kalman smoothing systems, shows that the backward algo- 
rithm can behave stably even when the tridiagonal system has 
null singular values. 

For y g R nxn we use the notation |v| for the operator norm 
of the matrix v; i.e., 



|v| sup{|vw| : w G R n , \w\ = 1} 
Theorem 5.2: Suppose that the matrices c^ and are given 



by 



c k = -q k a k 



where each G R nxn is positive definite, each e R nxn is 
positive semi-definite and a^ + \ = 0. It follows that <4 — g^T 1 is 
positive semi-definite for all k. Furthermore, if a is a bound 
for \qk\, \q^\, \uk\, and |%|, Then the condition number of dk 
is bounded by a 2 + a 6 . 

Proof: We note that = q~^ 1 + un so this conditions 
bound holds for k = N. Furthermore d^ — q^ 1 = so positive 
semi-definite assertion holds for k = N. 

We now complete the proof by induction; i.e., suppose 
<4 + i —q~/~+i is positive semi-definite 



dk 



h-cl+xd^Ck+i 



= u k + a M q M a k +i -a M q M d M q M a M 

= "^ + ^+1^+1 [ a k+i ^+1^+1 

The assumption that <4 + i — is positive semi-definite 
implies that that qk+\ — is positive semi-definite. It now 
follows that dk — q^ is the sum of positive semi-definite 
matrices and hence is positive semi-definite which completes 
the induction and hence proves that d^ — q^ 1 is positive semi- 
definite. 

We now complete the proof by should the condition number 
bound holds for index k. Using the last equation above, we 
have 



\d k 



< 



mi 



l4+il 2 k*+il 2 k*+il 



< a + a 



Hence the maximum eigenvalue of dk is less than or equal 
a + a 5 . In addition, since d^ — q^ 1 is positive semi-definite, 
the minimum eigen-value of dk is greater than or equal the 
minimum eigen-value of q^ 1 , which is equal to the reciprocal 
of the maximum eigen-value of qk- Thus the minimum eigen 
value of dk is greater than or equal 1/a. Thus the condition 
number of dk is bounded by a 2 + a 6 which completes the 
proof. ■ 

VI. Numerical Example 

To put all of these ideas in perspective, we consider a toy 
numerical example. Let n = 1 and N = 3, let q% = 1 for all k, 
and let <2k = 120. Then the system a T qa in ( |II.1| ) is given by 

"14401 120 ' 
120 14401 120 , 
120 1 

and its minimum eigenvalue is 4.8 x 10 -9 . To understand what 
goes wrong, we first note that the minimum and maximum 
eigenvalues of all the cik's except the last one coincide in 
this case, so the general condition of Corollary |2.1| will hold 
everywhere except at the last coordinate. 

Now that we suspect the last coordinate, we can check the 
eigenvector corresponding to the minimum eigenvalue: 

v min = [0.001 -.008 1] T . 

Indeed, the component of the eigenvector with the largest norm 
occurs precisely in the last component, so we are in the case 
described by Corollary |2.1| In order to control the condition 
number of this system, one option is to make the (3,2) and 



(2,3) coordinates less than 1 in absolute value. This changes 
the original system, but only at the last time point. Suppose 
we decide this is acceptable, and take these values to be 0.9 
instead of 120. The new system has lowest eigenvalue 1. If we 
expand the system (by increasing the size of {ak}, we 
find this stabilization technique works regardless of the size, 
since the last component is always the weakest link when all 
ak are equal. 

Next, suppose we apply the forward Thomas algorithm to 
the original (badly conditioned) system. We would get blocks 

d x = 14401, d 2 = 14400, d 3 = 4.8222 x 10~ 9 . 



This toy example shows Theorem 3.2 in action — indeed, the 



eigenvalues of the blocks are bounded above and below by 
the eigenvalues of the full system. Unfortunately, in this case 
this leads to a terrible condition number, since we are working 
with an ill-conditioned system to begin with. 

Now, suppose we apply the backward Thomas algorithm. 
We will get the following blocks: 

^3 = 1, di = 1 • d\ = 1 . 

Note that even though the blocks are stable, this does not mean 
that the backward algorithm can accurately solve the system 
a T qax = b. In practice, it may perform better or worse than 
the forward algorithm for ill-conditioned systems. 

Nonetheless, these results show the backward block Thomas 
has an important advantage: if it must be applied many times 
in an iterative context, it may happen that systems are badly 
conditioned at initial iterations, but the conditioning improves 
later on. In this context, the backward block Thomas it will 
run where the forward block Thomas will not. This was 
exactly the case in (5), where the introduction of the backward 
block Thomas stabilized the initialization procedure; for later 
iterations, the measurement terms came in for later iterations 
to improve the conditioning number of the system. 

VII. Conclusions 

In this paper, we have characterized the condition number 
of tridiagonal systems that arise in Kalman smoothing, see 



theorems 2.1 and 2.2 The analysis revealed that last block 
of the system can have a very strong effect on the system 
overall. In the Kalman smoothing context, it is important to 
note that the numerical stability conditions in theorems 2.1 
and |2.2| do not require ak in pi. 1| ) to be invertible. In fact, 
they may be singular, which means the process matrices G^, 
or derivatives of nonlinear process functions g£\ do not have 
to be invertible. 

We then showed that any well-conditioned symmetric block 
tridiagonal system can be solved in a numerically stable 
manner by the forward block Thomas algorithm. Applied to 
the Kalman smoothing context, we explicitly showed that the 
forward block Thomas algorithm is equivalent to the RTS 
smoother. Finally, we introduced a new method, the backward 
block Thomas algorithm, which shares the numerical stability 
properties of the forward block Thomas algorithm, and can 
remain stable even when applied to ill-conditioned systems. 



This property makes the backward Thomas algorithm prefer- 
able in an iterative context, where on occasion one may run 
into badly conditioned systems. The computational effort of 
the backward algorithm is equivalent to that of the forward; 
both scale linearly with the size of the tridiagonal system. 

Taken together, these results provide insight into both 
numerical stability and algorithms for Kalman smoothing 
systems. Block tridiagonal systems arise both in the classic 
setting, and also in many new methods, that can handle 
nonlinear process/measurement models, robust penalties, and 
inequality constraints. 
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